Numerical implementation of some reweighted path integral methods 
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The reweighted random series techniques provide finite-dimensional approximations to the quan- 
tum density matrix of a physical system that have fast asymptotic convergence. We study two 
special reweighted techniques that are based upon the Levy-Ciesielski and Wiener-Fourier series, 
respectively. In agreement with the theoretical predictions, we demonstrate by numerical exam- 
ples that the asymptotic convergence of the two reweighted methods is cubic for smooth enough 
potentials. For each reweighted technique, we propose some minimalist quadrature techniques for 
the computation of the path averages. These quadrature techniques are designed to preserve the 
asymptotic convergence of the original methods. 
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I. INTRODUCTION 

Random series representations are among the most 
general starting points for the derivation of finite approx- 
imations to the Feynman-Kag formula^ The Feynman- 
Kag formula computes the density matrix of a quantum 
system as an infinite dimensional stochastic integral 2,3 
and constitutes the main ingredient for the computation 
of quantal thermodynamic properties by Monte Carlo 
path integral techniques^ Current research is focused on 
the development of rapidly convergent finite-dimensional 
approximations to the Feynman-Kag formula. In decid- 
ing upon the efficiency of a certain method, one has to 
take into account both the number of path variables re- 
quired and the computational time necessary to evaluate 
the various quantities involved. 

One special group of path integral methods, which we 
shall refer to as direct methods, have the property that 
they require knowledge of the intermolecular potential 
only for the computation of the density matrix or the 
partition function of a physical system. Numerical al- 
gorithms based upon such methods have the advantage 
that they can be used as "black boxes." Many times, the 
potential energy surface of a system is a rather compli- 
cated function of the physical coordinates and there are 
no computational alternatives to such direct methods. It 
is therefore important to look for those direct methods 
that have the fastest asymptotic convergence and which 
minimize the number of calls to the potential. 

The most widely used direct method is the trapezoidal 
Trotter discrete path integral (DPI) mcthodi^ It is usu- 
ally derived by means of the Trotter product rule and 
an appropriate short-time high-temperature approxima- 
tion. The formal asymptotic convergence of the trape- 
zoidal Trotter DPI method and of related DPI techniques 
was extensively studied by Suzuki^ and was found to 
be 0(l/n 2 )£ As shown in Ref. such discrete path 

integral techniques can also be interpreted as direct dis- 
cretizations of the Feynman-Kag formula by means of 
certain quadrature rules. As a consequence, the dimen- 



sionality of the integral approximations for the density 
matrix is equal to the number of quadrature points. 

A second approach to constructing finite-dimensional 
approximations for density matrices is based on the ran- 
dom series representations of the Feynman-Kag formula. 
Special modifications of these random series by means 
of the so-called reweighted techniques have been initially 
considered on an intuitive basis in Refs. and De- 
pending upon the specific series, these modifications have 
been shown numerically to have asymptotic convergence 
equal to or faster than that of the trapezoidal Trotter DPI 
method. More recently, the correct definition and the 
convergence properties of the reweighted techniques have 
been shown to depend upon the asymptotic convergence 
of the partial averaging method>ii In fact, once the con- 
vergence of the partial averaging method has been fully 
understood?^ it became possible to develop reweighted 
techniques that have optimal asymptotic convergence^ 
The numerical implementation of two optimal reweighed 
methods that are based upon the Levy-Ciesielski and the 
Wiener-Fourier series, respectively will be thoroughly an- 
alyzed in the present article. 

Though their mathematical justification is rather in- 
volved, these two reweighted techniques have the advan- 
tages that they are direct methods and that they have cu- 
bic convergence for smooth enough potentials. However, 
as with the full Feynman-Kag formula, the path averages 
of the potential, which involve one-dimensional integrals 
on a compact interval, are assumed to be computed ex- 
actly. In practice, one is forced to make recourse to a 
certain quadrature scheme. Let us remember that for the 
trapezoidal Trotter DPI method the number of quadra- 
ture points is equal to the number of path variables. It 
follows that in order to preserve the overall cubic conver- 
gence of the reweighted methods, one needs to utilize a 
quadrature technique which uses a number of quadra- 
ture points proportional to the number of path vari- 
ables and which preserves the asymptotic convergence 
of the reweighted techniques. It is within the scope of 
the present article to show that such quadrature schemes 
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exist and to construct minimalist examples for either of 
the two reweighted techniques presented. 

The main tasks of the present article are a) to re- 
view in a transparent fashion the construction of the 
reweighted techniques, b) to present numerical exam- 
ples demonstrating the cubic convergence of the meth- 
ods, and c) to devise quadrature schemes that preserve 
the asymptotic convergence of the methods, while em- 
ploying a number of quadrature points proportional to 
the number of path variables. Based upon the numerical 
evidence presented as well as on the formal convergence 
properties of the methods, we eventually conclude that 
the reweighted techniques are the most efficient direct 
methods currently available. As such, these techniques 
constitute the methods of choice for use in path integral 
simulations. 



II. REWEIGHTED RANDOM SERIES 
METHODS 

In this section, we give a short review of the path in- 
tegral methods we study in this paper. According to 
Ref. (Q), the most general random series representation 
for the Feynman-Kag formula can be constructed as fol- 
lows. Assume given {Xk(r)}k>i a system of functions on 
the interval [0, 1], which together with the constant func- 
tion Ao(t) = 1 makes up an orthonormal basis in L 2 [0, 1]. 
Let 

Afe(t) = / X k (u)du. 



The primitive methods of order n are obtained by trun- 
cating the series appearing in the above formula to the 
rank n: 



p£ r (s,a';/?) 

p fp (x,x';(3) 
x exp 



d/i(ai 



dp,(a n ) 



(3 I V 

o 



x r (u) + a } j a k Afc (u) 



du 



k=l 



(4) 



It is known that the primitive methods can achieve 
at most 0(1/ n) asymptotic convergence^ and therefore 
they are of limited use. This is true especially because 
there are alternative methods capable of attaining faster 
asymptotic convergence without any additional work. 
Such methods include the so-called reweighted random 
series techniques, which will be discussed in the remain- 
der of the present section. 

According to Ref. lfTl|) . a reweighted method con- 
structed from the random scries YlkLx a k^k(u) is any 
sequence of approximations to the density matrix of the 
form 



Pf p (x,x';P) 



xexp<^ -p V 



d/i(ai 



d^(a qn+p ) 

qn+p 



c r (u) + a a k^n,k{u) du>, (5) 



k=l 



where q and p are some fixed integers, where 



A„,fc(u) = Afc(u) if 1 < k < n, 



(6) 



If fl is the space of infinite sequences a = (ax, et2, . . .) and and where 



dP[a] = n M 



Ok ) 



(1) 



k=l 



is the probability measure on ft such that the coordi- 
nate maps a — > a/c are independent identically distributed 
(i.i.d.) variables with distribution probability 



1 



/2tt 



- Q ?/ 2 da ; ; 



(2) 



then the Feynman-Kag formula has the form 
p(x,x';p) r ( 



pf p (x,x';/3) 7fi dP ^ eXP 

OO 

+ (7^tttA|-(u) 

fc=l 



P I V 

o 



x r (u) 



du 



(3) 



In the above equation, x r (u) = x + (x' — x)u is the ref- 
erence path, a = (h 2 ft /mo) 1 / 2 , and pf p (x,x'; (3) denotes 
the density matrix for a similar free particle. For a mul- 
tidimensional system, the Feynman-Kag formula is ob- 
tained by employing an independent random series for 
each additional degree of freedom. 



qn+p 

E 

k=n+l 



A„, fc (u) 



k=n+l 



(J) 



Notice that the expression given by Eq. (J5J is quite 
general and comprises the primitive methods by setting 
A„,fc(u) = for k > n + 1. Of course, in this case Eq. Q 
no longer holds. 

The reweighed techniques are not uniquely defined by 
the series representations on which they are based. There 
are an infinite number of reweighted techniques associ- 
ated with a given series. Some arc better, others arc 
worse. However, they share the common property that 
their asymptotic convergence is known»li The respective 
convergence theorems constitute definite criteria for the 
selection of optimal reweighted techniques. Two such op- 
timal methods will be presented in the remainder of the 
present section. 



A. A Wiener-Fourier reweighted technique 

The Wiener-Fourier series representation is based on 
the system of functions {A/c(r) = V2 cos(fc7rr); k > 1}, 
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which together with the constant function makes up a 
complete orthonormal system. In this case, we have 



A fe (i) = J \/2cos(fc7rT)dr = -y/J 



2 sin(fc7rt) 



k 



It has been showni that the Wiener-Fourier series rep- 
resentation is the optimal scries representation as to the 
minimization of the number of variables used to param- 
eterize the paths. This has been argued to be true of the 
primitive and partial averaging techniques, but it is also 
true of the reweighted techniques. 

A reweighted method based upon the Wiener-Fourier 
series representation and having superior asymptotic con- 
vergence can be constructed as follows iii We let q = 4 
and p — in Eq. JSJ) and define 



r n(u) 



l j n (u, U) 

7°(u, u)' 



(8) 



where 



2 sm(kiru) 2 
7 n (u,w) = -2 2^ 



fc=n+l 



2 ^ sin(/c7ra) 2 

= u(i - to - -j 2^ ■ 



2a n 



k=l 



in 



k 2 



k=n+l 



(9) 



(10) 



and 



Xn = L ^ u ^ du = V 2 £ Y 2 = \-Y^h- 

J " k=n+l fc=l 



Next, we let 



A„, fe (u) - Vi , 



2 sin(fc7ru) 



(11) 



(12) 



for 1 < k < n and 



2a n 



A„. fe (u) = \ ——r n (u) sin(fc7ru) (13) 
V on 

for n < k < An. One easily verifies that 



in 



£ k n.k(u) 2 = ^r„(u) 2 £ sin(fc 



3n 



7TU 



= 7„(u,u)= 2J A fc (w) 2 . 



Clearly the method defined by Eqs. JSJ, 112|) . and 
(|13|1 is a reweighted technique derived from the Wiener- 



Fourier series representation. Its explicit formula is 

p™ F (x,x>;f3) 



p f p(x,x';f3) 
x exp • 



= / d/i(ai). 

x r (u) + q-y^afc 



/? / y 

o 



d^(a 4 „) 

2 sin(fc7ru) 



fc=i 



2a, 



-o-A/ -5— r, 
3n 



,(m) 2J a/cSin(fc7ru) 

fc=n+l 



du 



(14) 



The reweighted technique defined above will be de- 
noted by RW-WFPI even if, strictly speaking, this 
acronym should denote the class of all Wiener-Fourier 
reweighted methods. Even if it might not be the best 
reweighted technique in its class, this particular RW- 
WFPI technique has good asymptotic convergence, as 
described by the following convergence results. Assume 
that V(x) has finite Gaussian transform and has sec- 
ond order Sobolev derivatives [for the exact conditions 
in which the results hold, the reader should consult 
Ref. (11)]. Then, 



lim n 3 [p(x,x';0)-p^ F (x,x , ;(3)] 



h 2 p 3 



NQ7T 4 mo 



where 



xp(x,x';p) \V'(x) 2 + V'(x') 2 ], 



N Q = 15.0045 . 



(15) 



(16) 



In the case of a d-dimensional system, the convergence 
constant becomes 



lim n 3 [p(x,x';/3)-p™ F (x,x';p)] = 



h 2 (3 3 



<p(x,x';/3) I £ 



IdtVjx)} 2 + [djVjx')]* 
m i 



(17) 



where diV (x) denotes the first order partial derivative 
dV{x)/dx t . 



B. A Levy-Ciesielski reweighted method 

Another important representation of the Feynman- 
Kag formula is based upon the so-called Levy-Ciesielski 
seriesjiS which is constructed as follows. For k = 1, 2, . . . 
and j — 1,2,..., 2 fc ~ 1 , the Haar function j is defined 
by 

r 2 (*-i)/a ) te[{l-i)/2 k ,l/2 k ] 

fk,j(t) = I -2W, t e [l/2 k , (I + l)/2 k ] (18) 
[ 0, elsewhere, 

where I = 2j — 1. Together with fo = 1, these func- 
tions make up a complete orthonormal basis in L 2 ([0, 1]). 
Their primitives 
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2 (*-i)/a[ t — (Z — l)/2*], te [(l-l)/2 k ,l/2 k ] 
Fkd(t) = { 2( fc - 1 )/ 2 [(; + l)/2 fe - t], te [l/2 k , {I + l)/2 k ] (19) 
0, elsewhere 



1.0 




u 



FIG. 1: A plot of the renormalized Schauder functions for 
the layers k — 1, 2, and 3 showing the pyramidal structure. 

are called the Schauder functions. The Schauder func- 
tions resemble some "little tents" that can be obtained 
one from the other by dilatations and translations. More 
precisely, we have 

F kA (u) = 2-( fc - 1 )/ 2 F 14 (2 fe - 1 U ) (20) 

for k > 1 and 

F kij (u)=F kil (u- 3 -^) (21) 

for k > 1 and 1 < j < 2 fe ~ 1 . The scaling relations above 
have to do with the fact that the original Haar wavelet 
basis makes up a multiresolution analysis of L 2 ([0, 1]) or- 
ganized in "layers" indexed by k. If we disregard the 
factor 2' fe ~ 1 )/ 2 j the Schauder functions make up a pyra- 
midal structure as shown in Fig.^ Let us notice that any 
reasonable Levy-Ciesielski primitive method should use 
a total of n = 2 k — 1 Schauder functions, corresponding 
to k complete layers. We shall assume that n = 2 k — 1 
for the remainder of the section. 

According to Ref. |TlT) . a Levy-Ciesielski reweighted 
method having superior asymptotic convergence can be 
constructed with the help of 3 ■ 2 k = in + 3 additional 
functions organized in two layers of 2 k and 2 k+1 func- 
tions, respectively. These additional functions are con- 
structed as follows. First, we define a function r(u) which 
is zero outside the interval (0, 1) and equal to 

r u(i-u) i 1/2 

\ w 2 F hl {u) 2 + (4 - 2u- 2 ) [F 2jl ( M ) 2 + F 2 , 2 (u) 2 ] J 

(22) 



C (u) 




FIG. 2: Shapes of the functions Cq(u), Lq(u), and Ro(u) 
utilized for the the construction of the reweighted technique. 

on the interval (0,1). The constant w appearing in 
Eq. J23 has the value 

w = 0.62258 .. . (23) 

Next, we define the functions: 

C {u) = wr(u)Fi,i(u), 



L (u) = y/4-2w 2 r(u)F 2>1 (u), 

and 

R (u) = ^/A-2w 2 r{u)F 2 ^{u). 

The functions Co(tt), Lq(u), and Ro(u) are plotted in 
Fig. □ They vanish outside the intervals [0, 1], [0, 1/2], 
and [1/2, 1], respectively. From them, by dilatations and 
translations, we construct 

C M (u) = 2- k ' 2 C (2 k u) , L kA (u) = 2- k ' 2 L (2 k u) , 

R ktl {u) = 2- k / 2 R (2 k u), (24) 

as well as 

C ktj (u)=C ktl [u~(j~l)/2 k ], 

L k<j (u) = L kil [ u -(j~l)/2 k ] , (25) 

R k ,j(u) = Rk,i [u - (j - l)/2 k ] , 

for l<j< 2 k . 



5 



In order to have a more compact notation, we set (re- 
member that n = 2 k — 1) 

F l ^(u) = F l , j (u) 

for 1 < I < k and 1 < j < 2 i_1 , 



p f p(x,x';f3) 



where [2 /_1 u] is the integer part of 2 l ~ 1 u. 

The particular reweighted Levy-Ciesielski method de- 
scribed in the present subsection will be denoted by RW- 
LCPI despite the fact that this acronym should refer 
to the whole class of reweighted techniques based upon 
the Levy-Ciesielski series. The method was proven to 
have o(l/n 2 ) convergence for potentials having first order 
Sobolev derivatives. It was conjectured that the asymp- 
totic rate of convergence may reach (9(l/n 3 ) for poten- 
tials having second order Sobolev derivatives. As a nu- 
merical advantage, we notice that the computation of the 
path at a point u involves only k + 2 = 2 + log 2 (n + 1) 
multiplications and k + 1 = 1 + log 2 (n + 1) additions. 



for 1 < j < 2 k , and 

rW f „\ _ / L ktP {u), if j = 2p- 1, 
W J -1 Rk, P (u), if J = 2p, 

for 1 < j < 2 fe+1 . Then, the reweighted technique has 
the explicit form 



(26) 



where po (x,x'; (3) stands for a so-called short-time high- 
temperature approximation. 



In this paper, we shall employ for comparison purposes 
the trapezoidal Trotter-Suzuki short-time approximation 
given by the formula 



j d« M ...J da fc+2 , 2 * +1 {2,y^ + ^ exp ( - \ £ £ 



2 ^— ' ^— ' ~ Zj 

'=1 3=1 



x exp < — j V 



fe+2 



I 



du 



III. DISCRETE PATH INTEGRAL METHODS 



Usually, by a discrete path integral method one under- 
stands an approximation to the density matrix having 
the Trotter product form 



p? FI (x,x';f3) = / dx 1 . . . [ dx n p Q (x,xi; - 

Jr Jr V n + 1 



Po x n ,x 



n + 1 



Po (x,x';/3) = pf p (x > x';f3)exp 



V(x) + V(x') 



(28) 

The resulting method will be denoted by the acronym 
TT-DPI. It has been shownifithat for n = 2 fe -l, the TT- 
(27) DPI method admits the following fast implementation 



p fp (x,x';(3) 



= / dai, 



da 



fe,2 fc -! 



(2tt 



-n/2 



1=1 i=l 



k T 



(29) 



exp < -P^uiiV 



i=0 



.(Uj) + 0-^ ^,[2'- 1 u j ]+l('"i) a J,[2 l - 1 « 4 ]+l 



2=1 



where it, = 2 fc i for < i < 2 h and 

2 -(fe+i) ! if i e {0,2 fe }, 
< 2- fe , if 1< i < 2 k - 1. 



The advantage of this scheme is the log(n + 1) scaling of 
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the number of operations one needs to perform in order 
to compute the path at a discretization point Uj. This is 
accomplished without the use of a fast Fourier transform 
routine. The asymptotic convergence of the method is 
known to be 0(l/n 2 ) for smooth enough potentials. 6 

A second discrete path integral method we consider in 
this paper is related to the RW-LCPI method. It has 
been shown that for n = 2 k — 1, the RW-LCPI method 
can be put in the Trotter product formii 

Pn C {x,x';f3) = J dx!...J dx n p^ c (x, x x ; 



L( I f @ 

■ Po I X n i X i 



where 

p% c (x,x';(3) _ 1 



p fp {x,x';l3) (2tt) 3/2 



n + 1 J ' 



-{al+a 2 2 +at)/2 



(30) 



xexp|-/3y V[x + (x' - x)u + aiaC (u) (31) 
+a2<jLo(u) + a^aRo(u)]du >daida2da3. 



However, for arbitrary n, the right-hand side of Eq. H30|) 
defines a new quantity which wc shall denote by 
Pn°{x, x'; 0), too. Because the n = 2 k — 1 subsequence 
has o(l/n 2 ) convergence, it is not farfetched to assume 
that the entire p\ (x, x';0) sequence has o(l/n 2 ) con- 
vergence. We shall verify this assumption later with the 
help of some numerical examples. 

One of the main advantages of the discrete path inte- 
gral methods consists of the fact that for low dimensional 
systems the evaluation of the density matrix and related 
properties can be performed accurately by means of the 
numerical matrix multiplication (NMM) methodpi^i^ We 
shall use the NMM method to compute n-th order ap- 



proximations to the partition function of the type 



for one-dimensional systems. The main steps of the al- 
gorithm are as follows. First, one restricts the system to 
an interval [a, b] and considers a division of the interval 
of the type 

Xi = a + i(b - a) /M, < i < M. 

Next, one computes and stores the symmetric square ma- 
trix of entries 



b — a 
M 



Po° I x h x ji 







71 + 1 



, 0<i,j<M. 



The value of the partition function can then be recovered 
as 

By computer experimentation, the interval [a, b] and the 
size M of the division are chosen such that the compu- 
tation of the partition function is performed with the 
required accuracy. A fast computation of the powers 
of the matrix A can be achieved by exploiting the rule 
j^m+n _ (^4»«)»\ p or more details, the reader is referred 
to the cited literature. 

The use of the NMM technique together with the 
trapezoidal Trotter approximation raises no special com- 
putational problems. However, for the DPI method 
derived from the Levy-Ciesielski reweighted technique, 
there is the problem of evaluating the short-time approx- 
imation given by Eq. (|31l) . Noticing that the function 
Lq{u) vanishes outside the interval [0, 1/2] while the func- 
tion Rq{u) vanishes outside the interval [1/2, 1], one may 
recast Eq. <|31ll in the form 



p fp (x,x';(3) 



d/x(ai) 



1/2 



dp,{a%) exp i — (5 J V[x + (x! — x)u + ai<rCo(u) + a2<jLo(u)]du 

/ d/i(c(2) exp < — (3 I V[x + (x' — x)u + aicrCo(u) + a,2crRo(u)]du 
Jr I Ji/2 



(32) 



which is computationally less expensive because it no 
longer involves a three dimensional Gaussian integral. 
The Gaussian integrals appearing in Eq. I|32[l [remember 
Eq. @] can be evaluated by means of the Gauss-Hermite 
quadrature techniqueii^ For the purpose of establishing 
the asymptotic convergence of the partition functions, 
we have found that a number of 10 quadrature points for 
each dimension is sufficient. This is so because the er- 



rors due to the Gauss-Hermite quadrature approximation 
quickly vanish as (3/(n + 1) — > 0. 
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IV. PATH INTEGRAL QUADRATURE 
TECHNIQUES 

The basic assumption utilized in the construction of 
the reweighted techniques is that the path averages of 
the type 



V 



x r (u) 



qn+p 

a 2J a k A nt k(u) 
fc=i 



du 



are computed exactly. In practice, however, one has to 
perform the integrals approximately by quadrature tech- 
niques. In this section, we shall give a positive answer to 
the question: Are there any quadrature rules that pre- 
serve the fast asymptotic convergence of the reweighted 
methods while employing at most q'n + p' quadrature 
points for some fixed integers q' and p'l 

It is not within the scope of the present paper to solve 
this problem in a mathematically rigorous way. Rather, 
we shall use numerical examples to demonstrate our find- 
ings. The key observation toward the solution of the 
quadrature problem is the fact that the asymptotic rates 
of convergence of the reweighted and partial averaging 
techniques depend solely on the first and the second order 
derivatives of the potential, through their square. Exam- 
ples are furnished by the convergence constants of the 
partial averaging Wiener-Fourier method (Theorem 4 of 
Ref. Il2f) . the partial averaging Levy-Ciesielski method 
(Theorem 4 Ref. OH, or the Wiener-Fourier reweighted 
method discussed in Section II. A. 

This observation suggests that any quadrature rule ca- 
pable of integrating exactly or to a good approximation 
any quadratic potential will preserve the cubic conver- 
gence of the respective methods for all smooth enough 
potentials (potentials having continuous second order 
derivatives). In the present section we shall focus our 
attention on finding the fastest quadrature techniques 
that integrate the quadratic potential. The fact that the 
cubic convergence of the reweighted methods extends to 
all sufficiently smooth potentials is then demonstrated in 
the following sections with the help of some numerical 
examples. 



Minimalist quadrature rule for the RW-LCPI 
method 



In this subsection, we shall find an explicit quadrature 
technique that integrates exactly all the quadratic poten- 
tials for the RW-LCPI method and its discrete path inte- 
gral version. We look for a quadrature technique that uti- 
lizes a minimal number of quadrature points and which 
also preserves the intrinsic symmetries of the system of 
functions utilized for path expansion. 

We begin by noticing that on each interval 
[(j - l)2- fc , j2" fc ] the paths appearing in the RW-LCPI 
formulation are linear combinations of the functions 1, u. 



C n ,j{u), L n j(u), and R n j(u). That is, they are of the 
form 

c\ + c 2 u + czC n> j(u) + CiL n j(u) + c 5 R, hj (u) 

with some appropriate values for the constants c\, . . . , c&. 
This is so because the functions Fi jP are linear on the 
intervals [{j - l)2~ fc , j2~ k ] for all 1 < I < k and 1 < p < 
2 l ~ 1 . Using this observation and the decomposition 



V 



k+2 



du 



k+2 

a ^,[2 i - 1 u] + l 



E / v 



i=i 



du 



together with Eqs. (f2~l)l and f2*K|l . it is not difficult to 
establish the following proposition. 

Proposition 1 Assume that the quadrature rule speci- 
fied by the the points — uq < u\ . . . < u q = 1 and the 
nonnegative weights wq, . . . , w q integrates exactly all the 
expressions of the form 



fv[ Cl 

Jo 



+ c 2 u + c 3 C (u) + caL q {u) + c 5 R (u)] du 



for a given class of continuous functions V(x). Then the 
quadrature rule specified by the q2 k + 1 points 



u' p = (u p - q \p/ q] + [p/q]) 2- 
and the corresponding weights 



forO<p< q2 k 




n — k 
W p-q[p/q] Z 



for p = 

forp = q2 k , 

forp = qj, 0<j <2 k , 

otherwise 



integrates exactly all the expressions of the form 

k+2 

x r (u) + a^2a lA2 ,-i u]+1 F^l_ lu]+1 (u) 
1=1 



V 



du 



for the same class of potentials. 

While the straightforward proof of the proposition is left 
for the reader, we mention that an alternative demon- 
stration can be achieved by means of Eqs. i|30|l and 13111 . 
Again, [p/q] denotes the integer part of p/q. 

Proposition Q reduces the quadrature problem to the 
problem of integrating expressions of the type given by 
Eq. (|33|l . which are simpler. A further reduction can 
be achieved by noticing that the function Cq(u) is sym- 
metrical under the transformation u — > 1 — u, while the 
functions Rq(u) and Lq(u) transform one into each other. 
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Moreover, any linear combination cq + c\u remains a lin- 
ear combination of the functions 1 and u under the same 
transformation. These arguments allow us to restrict our 
search to the interval [0, 1/2] only, because the rest of the 
quadrature points on the whole interval [0, 1] can be ob- 
tained by symmetry. To integrate all the expressions of 
the form 

,1/2 

/ V [ci + c 2 u + c 3 C (u) + c 4 L (m)] du 
Jo 

(the function Rq(u) disappeared because it is zero on the 
interval [0, 1/2]) for all quadratic potentials, it is enough 
that the following system of equations be satisfied: 



Ei=0 



,1/2 



,,-n — /( 
Ei=0 W i U i 

E4 2 
i=0 W i U i 

Ei=0 w i C o(Ui) ■■ 

Ei=0 WiC (Ui) 2 
,4 



Jo 



1/2 



udu, 



r 1 / 2 , 
Jo 

) = 

) 2 



' u 2 du, 

ti /2 C (u)du, 
= J 1/2 Co(ufdu, 

"1/2, 



(34) 



Ei=o w i L o( u i) = Jq /2 L (u)du, 
Y?i=o WiUiL^m) = fy /2 uL (u)du, 



Notice that the condition 



A f 1/2 
2_^WiLo{ui) 2 — / L (u) 
Jo 



'du 



(35) 



was not introduced in the above system of equations be- 
cause it is not independent of the others. More precisely, 
we have the relation 

C a (u) 2 + L {u) 2 =u(l-u), 

which shows that the equality (|35|l is automatically sat- 
isfied provided that the first three and the fifth equalities 
from the above system hold. 

The system given by Eq. I)34|) contains 9 equations and 
this is why the minimal number of quadrature points for 
the interval [0,1/2] is 5. Consequently, the number of 
unknown variables is 10 and so we can arbitrarily fix 
one of the variables u,. We have chosen to set uq = 
0. With this choice, the system of equations described 
by Eq. I)34|) can be solved by means of the Levenberg- 
Marquardt algorithm as implemented in Mathcad. 16 The 
computed quadrature points and weights together with 
their extension to the whole interval [0, 1] are shown in 
Table HI Luckily enough, the weights are all positive and 
therefore the quadrature scheme is numerically robust. 



It is instructive to verify numerically the findings of the 
present subsection for a one-dimensional system made up 
of a particle of mass mo = 1 moving in the quadratic po- 
tential V(x) = m uj 2 x 2 /2 of frequency uj = 1. We use 
atomic units (therefore, h = 1) and we set (3 = 10. As 
discussed in the preceding section, it is more convenient 
to study the convergence of the discrete path integral 
method based on the Levy-Ciesielski reweighted tech- 
nique because this can be performed by numerical matrix 
multiplication and it coincides with the corresponding se- 
ries representation for n = 2 k — 1. 

By numerical matrix multiplication, we compute a se- 
quence of partition functions ^m+it/^) an d then we com- 
pute the ratios 



as well as the quantities 



y^ c = m 2 In 



1 + 



Here, Z(f3) denotes the exact partition function. Assum- 
ing that the asymptotic convergence of the odd sequence 
1 (/3) can be described by the equation 



f?LC 
ri 2m+ 



i££+i(/3) = l + 



clc 



-LC 



(2m + 1) Q 



(2m + 1) Q+1 
1 



(2m + 1) Q+1 



(36) 



it is not difficult to show that the slope of as a func- 
tion of m converges to the convergence exponent a as 
fast as o(l/m)i Moreover, once the exponent a is deter- 
mined, one may compute the quantities 

4 c = (2m + irm [i^ +1 (/3)-l] 



and show that c^ +1 



converges to clc as fast 



as o(l/m). A similar discussion is true of the even 
sequence B^miP)^ though for the case of the Wiener- 
Fourier reweighted technique, the constant c'^q appearing 
in Eq. Il.'itill has usually different values for the odd and the 
even subsequences. To avoid the appearance of certain 
oscillations in our plots as well as to maintain a unified 
exposition of the methods employed, we study only the 
odd subsequence but mention that the even subsequence 
was also studied and produced the same results. 

The short-time high-temperature approximation uti- 
lized in the NMM technique has the expression 
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TABLE I: Quadrature points and weights for the RW-LCPI method. 



i 





1 


2 


3 


4 


5 


6 


7 


8 


9 


Ui 


0.000000000 


0.081468437 


0.232499762 


0.36369423 


0.468289398 


0.531710602 


0.63630577 


0.767500238 


0.918531563 


1.000000000 


Wi 


0.014247669 


0.149656564 


0.131280659 


0.130290452 


0.074524656 


0.074524656 


0.130290452 


0.131280659 


0.149656564 


0.014247669 



p f p(x,x';(3) 



d//(ai) 



/ d^(a 2 )exp<^ — fi'S~]wiV[x + (x' — x)ui + aiaCo(ui) +a2a-Lo(ui)]> 
/ d^(a 2 )exp^ - fi'S^WiV[x + (x - x)ui + aiaC${ui) + a2aRo(ui)] 



(37) 





FIG. 3: The convergence order of the RW-LCPI method 
depends upon the nature of the quadrature scheme employed. 



FIG. 4: The extrapolated value of the convergence constant 
for the RW-LCPI method is c L c = 10.32. 



where the one-dimensional integrals appearing in 
Eq. (l.'S2t have been transformed into finite sums with the 
help of the quadrature points tabulated in Table [I] As 
discussed in the preceding section, the Gaussian integrals 
are computed by Gauss-Hermite quadrature in 10 points 
for each dimension. We have also computed a sequence 
a'^ C for which the 10 quadrature points from Table [I] 
were replaced by the same number of Gauss-Legendre 
quadrature points^ 



As shown in Fig. |3 the sequence a^ +1 



... . a n? con - 
cv'^j , which converges 



verges to 3 as opposed to a' m+1 
to 2. This demonstrates that both the number of quadra- 
ture points and the nature of the quadrature technique 
play an important role for the value of the convergence 
order. In fact, additional numerical examples not pre- 
sented here show that doubling the number of quadrature 
points for the Gauss-Legendre method does not improve 
the convergence order. We have also tested the trape- 
zoidal quadrature scheme and found again a convergence 
order of 2. Therefore, the quadrature scheme is impor- 
tant and it cannot be chosen arbitrarily. We therefore 
recommend the use of the quadrature scheme developed 



in the present section for practical applications. In the 
following section, we will demonstrate by numerical ex- 
amples that this scheme preserves the cubic convergence 
of the RW-LCPI method for all smooth enough poten- 
tials. 

In Fig. EI we plotted the sequence t^+i — c m°' tne 
limit of which is the convergence constant clc ~ 10.32. 
This value of the convergence constant will be compared 
in the following subsection with the one for the Wiener- 
Fourier reweighted technique. As discussed in Section II, 
the Wiener-Fourier reweighted technique is expected to 
have a smaller convergence constant. 



B. Minimalist quadrature rule for the RW-WFPI 
method 



Theoretically speaking, we may use a strategy similar 
to the one employed for the RW-LCPI method in order to 
develop a quadrature technique which integrates exactly 
all the quadratic potentials. However, for the RW-WFPI 
method, this approach is rather problematic because it 
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requires solving large dimensional non-linear systems of 
equations of the type shown by Eq. Also, because 

the products of the functions utilized for path expansion 
present a strong linear dependence even after all func- 
tional relations are taken into account, the resulting sys- 
tem of equations is rather ill-conditioned. 

The alternative approach advocated in the present sub- 
section is an empirical one. Let us remember that the 
functions utilized in the path expansion are 1, u, sin(/c7ru) 
for 1 < k < ri, and r n (u) sin(fc7ru) for n < k < An. It 
is known that lim n ^oo r n (u) = l4i It follows that these 
path expansion functions as well as their products are 
highly oscillatory functions which may be well approx- 
imated by the Legendre orthogonal polynomials on the 
interval [0, 1]. This might be so because the latter ones 
are also highly oscillatory. This suggests that the Gauss- 
Legendre quadrature scheme (which we have found does 
not perform well for the RW-LCPI method) may be a 
good quadrature scheme for the RW-WFPI method. 

Before moving on to numerical examples, let us de- 
fine what we mean by "good" quadrature schemes. We 
require that the quadrature scheme preserve the asymp- 
totic convergence of the RW-WFPI method. More pre- 
cisely, both the order of convergence and the convergence 
constant should be preserved. Let us define 



c™ F =(2m + l) 



Z(J3) 



(38) 



Again, we study the convergence of the odd subsequence, 
but mention that the even subsequence was also studied 
and produced similar results. From Eq. I|15|l . one com- 
putes the following theoretical value for the convergence 
constant 



c WF = Jim_ c, 
1 



WF 



2 .. v:S 



Z{fi) 



N Q ir 4 m 
p(x,x;P)V'{x) 2 dx re 



-0.684. (39) 



The numerical evaluation of the sequence c^ F by ma- 
trix diagonalization has been discussed in Appendix C 
of Ref. IjllT) . Employing the method developed there, we 
have computed numerical approximations to the terms 
c^ F using 6n, 7n, and 8n Gauss-Legendre quadrature 
points, respectively (remember that n = 2m +1). The 
results plotted in Fig. clearly demonstrate that the con- 
vergence order is 3 for all three cases. However, one no- 
tices that if the number of quadrature points is larger 
or equal to 7n, the results are virtually indistinguishable 
and the computed values converge to one and the same 
limit. This limit is automatically the theoretical conver- 
gence constant. 

We can reinforce the findings of the preceding para- 
graph by studying the convergence of the sequence 



m Ic, 



WF 



L WF 



i. If the sequence is convergent, then 



lim c, 

m — >oo 



.WF 



C WF- 



□ - 6n points 

O - In and 8n points 




-1 



c%= -0.684 



10 



30 
m 



50 



FIG. 5: The convergence constant for the RW-LCPI method 
depends upon the number of quadrature points in the Gauss- 
Legendre scheme. Using In — 7(2m + 1) quadrature points 
or more is enough to recover the correct theoretical constant. 



20 



o 

fe'g 10 

o 



□ - 6n points 
O - In and 8n points 




10 



30 
m 



50 



FIG. 6: This plot shows the qualitative difference for go- 
ing from 6n to 7n or more quadrature points. In the latter 
case, the theoretical convergence constant of the RW-WFPI 
method is preserved. 



On the other hand, if the asymptotic dependence of 
with m is linear with a non-zero slope, 



m [c, 
then 



WF 



„th 
L WF, 



lim 



WF 



7^ C WF 



The results plotted in Fig. clearly demonstrate that 
the asymptotic convergence of the RW-WFPI method 
(meaning both the convergence order and the conver- 
gence constant) is preserved provided that the number 
of Gauss-Legendre quadrature points is 7n or larger. 

We conclude this section by performing a comparison 
between the rates of convergence of the RW-LCPI and 
RW-WFPI methods. For a given n, the former method 
employs An + 3 Gaussian parameters for path expansion, 
while the latter method employs An. The two numbers 
can be thought of as equal for large n. However, the 



11 



modulus of the convergence constant for the RW-WFPI 
method is 0.684, almost 15 times smaller than the con- 
vergence constant for the RW-LCPI method, which is 
10.32. As order of magnitude, the difference between the 
two constants at low enough temperatures can be under- 
stood as follows. As shown by Eq. (|39|) . the dependence 
of the convergence constant with f3 at low temperature 
is of the form 

c$ F (/?) oc /3 3 . 

On the other hand, assuming that the convergence order 
for the RW-LCPI method is 3 for all smooth enough po- 
tentials, the short-time high-temperature approximation 
given by Eq. I|31|) must be of the form 

pl c (x,x';f3)=p(x,x';0) [l + 0(/3 4 )] 

for small enough (3. This is so because the convergence 
order of the discrete path integral method described by 
Eq. H30|) and the order of the short-time approximation 
given by Eq. 1|31|) are always related by the following 
Suzuki formula™ 

p^ c (x, x'; P) = p(x, x'- /?) [1 + 0(/? 4 /n 3 )] . 

This formula shows that 

and explains why the Wiener-Fourier reweighted tech- 
nique performs better than its Levy-Ciesielski version at 
low temperature (large (3). 

V. MORE NUMERICAL EXAMPLES 

In this section, we present numerical evidence support- 
ing the idea that the convergence of the reweighted tech- 
niques described in the previous section is indeed cubic 
for smooth enough potentials even if only the proposed 
minimal quadrature techniques are employed. Strictly 
speaking, we will be able to verify our findings only 
for the RW-LCPI method by numerical matrix multi- 
plication. For the RW-WFPI method, we shall compute 
the partition functions of the systems studied by Monte 
Carlo integration and show that they are better than the 
partition functions obtained by means of the RW-LCPI 
method. Clearly, this provides only indirect evidence for 
the asymptotic convergence of the RW-WFPI method. 
As discussed in Section IV. B of Ref. (1), it is very diffi- 
cult to compute the partition functions by Monte Carlo 
integration with accuracy good enough to determine un- 
equivocally the convergence order of a method. 

Besides the convergence order, we shall also verify that 
the quadrature technique developed in Section IV. A pre- 
serves the convergence constant of the RW-LCPI method. 
Since we do not have an explicit expression for it, we de- 
termine the convergence constant numerically by using 
a large number of Gauss-Legendre quadrature points to 



compute the "exact" short-time approximation given by 
Eq. (J22J. Thus, we compute the sequences aj^ c and c^ c 
defined in Section III. A with the help of the 10 quadra- 
ture points from Table [I] but also the sequences a'^ 
and c^ c with the help of 500 Gauss-Legendre quadrature 
points. In the preceding section, we suggested that no 
matter how large the number of Gauss-Legendre quadra- 
ture points is, the convergence order of the resulting 
method is at most 0(l/n 2 ). However, for a fixed range 
of values of n, one can obtain essentially exact results for 
the quantities a'^ c and c^ c by employing a large number 
of Gauss-Legendre quadrature points. Numerical exper- 
imentation shows that the optimal number is 500 for the 
range of n employed in the present study (n — 2m + 1 , 
m= 1,2,..., 60). 

A. Quartic potential 

The first example we study is the quartic potential 
V(x) = x 4 /2. We set H = 1, m = 1, and /? = 10. With 
the help of the numerical matrix multiplication technique 
described in Section III. A, we compute the quantities 

and 

The exact partition function Z{j3) has been evaluated by 
variational methods and has the value 

Z(/3) w 4.982570651 • HT 3 . 

To compute the partition functions ^m+id^) by numer- 
ical matrix multiplication, we have restricted the prob- 
lem to the interval [—4, 4]. A division of this interval in 
M = 1024 slices is sufficient for Z 2 i4_ 1 (/3) to reproduce 
the first 9 digits of the exact partition function. 

As previously discussed, the slope of as a function 
of m converges to the convergence exponent a. Once the 
exponent a is determined, one may obtain the conver- 
gence constant by analyzing the slope of 

c^ c = (2m + l) Q m[i££: +1 (/3)-l]. 

A similar study is performed for the sequences a„ c and 
c^ c , the slopes of which converge to the exact conver- 
gence exponent and convergence constant, respectively. 

As demonstrated by Figs. □ and El the RW-LCPI 
method based on the quadrature technique developed in 
Section IV. A and the exact RW-LCPI method have vir- 
tually identical convergence properties. Both methods 
have cubic convergence and have equal convergence con- 
stants. This is strong evidence supporting the hypothesis 
that if a quadrature technique preserves the asymptotic 
convergence of a reweighted method for the harmonic os- 
cillator, then it preserves the asymptotic convergence for 
all smooth enough potentials. 
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FIG. 7: This plot actually contains two superimposed graphs 
showing that the slopes of a'^° and a^ converge to 3. 

150-, 1 1 1 1 1 r 



tion functions decrease up to a point (n = 32) and then 
increase back to the exact partition function. This is con- 
sistent with the fact that the convergence constant of the 
RW-WFPI method is negative. 

TABLE II: Partition functions for the quartic potential. The 
partition functions for the RW-WFPI method were computed 
by Monte Carlo integration in 1 billion points. We employed 
7n Gauss-Legendre quadrature points to compute the path 
averages. 



n 








2 


10.235 • 10~ 3 


10.688 ■ 10" 3 


(7.899 ± 0.003) ■ 10" 3 


4 


6.489 ■ 10~ 3 


6.538 ■ nn 3 


(5.518 ±0.003) ■ 10~ 3 


8 


5.393 • 10~ 3 


5.364 ■ 10" 3 


(5.047 ± 0.002) ■ 10~ 3 


16 


5.089 • 10" 3 


5.062 ■ 10~ 3 


(4.983 ± 0.002) • 10~ 3 


32 


5.009 ■ 10~ 3 


4.996 ■ 10~ 3 


(4.979 ± 0.002) ■ 10~ 3 


64 


4.989 ■ 10~ 3 


4.985 ■ 10~ 3 


(4.982 ± 0.002) ■ 10~ 3 


CO 


4.983 ■ 10~ 3 



100 




FIG. 8: This plot contains two superimposed graphs demon- 
strating that the slopes of c'^° and cj^ are virtually indistin- 
guishable. Therefore, the slopes converge to the same limit 
as m — + oo. 



Even if we cannot perform a similar analysis for the 
RW-WFPI method, it is instructive to compare the val- 
ues of the partition function with those computed by the 
trapezoidal Trotter and RW-LCPI methods. In Table 1TT1 
we have listed the partition functions Zj^_ 1 (0), Z^ ± (P), 
and Z^ F ((3) for n — 2 k , k = 1,...,6. By computing 
ZjJ_ 1 (f3), Z^ 1 (0), and Z^ F ((3), we ensure fair compar- 
ison as to the number of variables used to parameter- 
ize the paths. The first two methods utilize a number 
of An — 1 path variables, while the last one utilizes An 
variables. One notices that even if the Levy-Ciesielski 
reweighted method has cubic convergence, the improve- 
ment over the trapezoidal Trotter is not that great for 
small numbers of path variables. This is due to the fact 
that the convergence constant of the RW-LCPI method 
is large, being proportional to /3 4 . The Wiener-Fourier 
reweighted method, which has a convergence constant 
proportional to (3 3 , produces much better values for a 
similar number of path variables. The RW-WFPI parti- 



B. He cage problem 

The physical system^ we consider in this subsection 
consists of a particle trapped on a line between two atoms 
separated by a distance L. The particle is assumed to 
interact with the fixed atoms through pairwise Lennard- 
Jones potentials. The resulting cage is described by the 
potential 



V{x) = Ae 



12 



L 



12 



L 



if < x < L and V{x) — +oo, otherwise. The parame- 
ters of the system are chosen to be those for the He atom. 

o 

We set mo = 4 amu, e/k B = 10.22 K, a — 2.556 A, and 

o 

L = 7.153 A. At T = 5.11 K, which is the temperature 
utilized in the present computations, the system is prac- 
tically in its ground state. Numerical experimentation 
shows that in order to compute the partition functions 
ZJm+iO-O by numerical matrix multiplication, it is safe 
to restrict the problem to the interval [0.2L,0.8L] and 
consider a division of this interval in M = 1024 equal 
slices. By employing large values of m, the exact parti- 
tion function was found to be 

Z{(3) fa Z\^_ x (fi) = 1.668294002 . . . 

Figs. |5] and ^] demonstrate that the convergence of 
the RW-LCPI method for the He cage problem is cubic. 
Moreover, there is virtually no difference between the 
values computed using the quadrature technique of Sec- 
tion IV. A and those computed using 500 Gauss-Legendre 
quadrature points. Similar to the case of the quartic po- 
tential, the partition functions Z^-iift) and Zj^_ 1 (f3) 
are very close one to each other for small n even if 
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FIG. 9: This plot actually contains two superimposed graphs 
showing that the slopes of a'^ c and aj^ c converge to 3. This 
demonstrates that the convergence order of the RW-LCPI 
method for the He cage problem is cubic. 



60 




m 



TABLE III: Partition functions for the He cage problem. The 
partition functions for the RW-WFPI method were computed 
by Monte Carlo integration in 1 billion points and using In 
Gauss-Legendre quadrature points to compute the path aver- 
ages. 



n 


^4 T n T -l(/3) 


Z^i(0) 




2 


1.9161 


1.8917 


1.8084 ± 0.0003 


4 


1.7597 


1.7545 


1.7050 ± 0.0003 


8 


1.6975 


1.6958 


1.6739 ± 0.0003 


16 


1.6766 


1.6753 


1.6686 ± 0.0003 


32 


1.6705 


1.6698 


1.6681 ± 0.0003 


64 


1.6688 


1.6685 


1.6683 ± 0.0003 


CO 


1.6683 



and 

/ p(x,x;/3)V"(x) 2 dx 

are finite, then the reweighted methods have cubic con- 
vergence. (Remember that the expectation values of the 
squares of the first and second order derivatives control 
the convergence constants of the reweighted techniques.) 
In the case of the Lennard-Jones potential, exp[— j3V(x)] 
has continuous derivatives of any order and the density 
matrix has an exponential decay near singularities (see 
also the Appendix). Consequently, the expectation val- 
ues of the first and second derivatives are finite. We 
believe the cubic convergence of the reweighted methods 
extends to all potentials satisfying these criteria. 

VI. CONCLUSIONS 



FIG. 10: This plot contains two superimposed graphs demon- 
strating that the slopes of c^ c and cj^ c for the He cage prob- 
lem are virtually indistinguishable. Therefore, the slopes con- 
verge to the same limit as m — > oo, which demonstrates that 
nothing is lost by employing the quadrature technique devel- 
oped in Section IV. A. 

the asymptotic convergence of the methods is different. 
On the other hand, as seen from Table IIIII the parti- 
tion functions computed with the help of the Wiener- 
Fourier reweighted technique are significantly better. 
The asymptotic convergence of Z^ w (f3) is eventually 
from below, in agreement with the theoretical prediction 
for the exact Wiener-Fourier reweighted technique. 

The He cage problem is interesting because the 
Lennard-Jones potential lies outside the class of poten- 
tials for which the reweighted methods were proven to 
have cubic convergence^ The example studied suggests 
that if exp[— (3V(x)\ has Sobolev derivatives of up to the 
second order for all (3 > and if the quantities 

/ p(x,x;j3)V'(x) 2 dx 



In this article, we have discussed the numerical im- 
plementation of two special reweighted techniques based 
upon the Levy-Ciesielski and the Wiener-Fourier series, 
respectively. For both methods, we have developed min- 
imalist quadrature techniques capable of preserving the 
asymptotic convergence of the respective methods, while 
employing a number of quadrature points proportional to 
the number of path variables. Moreover, we have demon- 
strated by numerical examples that the asymptotic con- 
vergence of the methods is cubic for smooth enough po- 
tentials. The numerical evidence is consistent with the 
mathematical analysis previously performed by one of 
us-ii 

The findings of the present paper allow us to con- 
clude that the reweighted techniques have an asymptotic 
convergence better than that of the trapezoidal Trotter 
DPI method or other related direct methods. More pre- 
cisely, the reweighted Levy-Ciesielski method is the most 
efficient method if the calculation of the paths is the 
dominant computational factor (theoretically, this will 
always be the case for a large enough number of path 
variables n). The RW-LCPI method has cubic conver- 
gence and requires a number of operations proportional 



14 



to nlog 2 (n) for the computation of the the paths. The 
Wiener-Fourier reweighted method is the most efficient 
direct method as to the minimization of the number of 
calls to the potential. This is so because the convergence 
constant of the RW-WFPI method is smaller than that 
of the RW-LCPI method, especially at low temperature. 
Even if the number of operations necessary to compute 
the paths is proportional to n 2 , the RW-WFPI technique 
will most likely be the most efficient method for practical 
applications. Calls to the potential routine are usually 
expensive, whereas the calculation of the paths takes full 



advantage of the vector floating point units available in 
modern processors. 
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APPENDIX A: POTENTIALS WITH STRONG 
POSITIVE SINGULARITIES 

In this section, we shall discuss the class of poten- 
tials that are bounded from below and have strong pos- 



itive singularities at origin of the form r~ q with q > 0. 
This class of potentials contains the Lennard- Jones po- 
tential. We shall first discuss the case of three dimen- 
sional potentials and then go back to the one-dimensional 
Lennard- Jones potential discussed in Section V.B. We fo- 
cus mainly on the three-dimensional examples because 
these are of practical importance. As such, we consider a 
particle of mass too = 1 moving in a spherical potential 
of the form 

V(r) = U(r)+r 2 /2, 

where U(r) is a bounded from below potential that be- 
haves at origin as r~ q and is otherwise an infinitely differ- 
entiable function decaying to zero at infinity. The term 
r 2 is added to confine the system about the origin so that 
the resulting Hamiltonian has discrete spectrum. 

The question we ask is about the decay of the density 
matrix at the singularity r = 0. For the sake of having 
an example, we set U(r) = r~ q . Clearly, the classical 
diagonal density matrix exp[— f3V(r)] has an exponential 
decay to zero at the singularity any time q > 0. As the 
quantum diffusion takes place (we interpret j3 as a time 
parameter) , the system starts to explore the singularities 
of the potential and the diagonal density matrix tends 
to spread out. The strong decay to zero of the classical 
density matrix near singularities is attenuated and, as 
we shall see, if the singularities are not very strong, the 
quantum density matrix may even fail to decay to zero. 
On physical grounds, we expect the effect to get stronger 
in the limit f3 — > oo , when the system is also brought into 
its ground state eigenfunction. Therefore, it is not far- 
fetched to assume that the decay of the density matrix at 
the singularity is controlled by the decay of the ground 
state eigenfunction. For this reason, we shall focus our at- 
tention on the behavior of the ground state eigenfunction 
near singularities. We remind the reader that the ground 
state eigenfunction satisfies the Schrodinger equation 



1 l__d 

2 r 2 dr 



r 2 — ^ (r) 
dr 



V(r)*(r) = E $(r) 



and that it is the only positive eigenfunction (one uses 
this property to easily check if an eigenfunction is the 
ground state eigenfunction). 
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It has been discussed in Section II. B of Ref. l(l2|) that if 
q < 2, then the potential V(r) is a Kato-class potential. 
In this case, it is known that the density matrix does not 
decay to zero at the singularity. A quick way to verify 
this is to notice that if q < 2, then the potential V(r) has 
finite Gaussian transform. In this conditions, Theorem 4 
of Ref. © says that 



3 „ TD>3 



V[x r (u) + aB°] < oo, V (a;, x') £ W x 



and for all j3 > 0. The last inequality together with 
Jensen's inequality implies 



-£|^ =E exp{-^V M «) + .Ba 



> exp |-/3E J V[x r {u) + aB^]j > 0, 

which proves our assertion since the density matrix of a 
free particle is strictly positive. 

On the other hand, if q = 2, the density matrix has 
polynomial decay to zero at r = 0, as shown by the decay 
of the (unnormalized) ground state eigenfunction 



*o(r) 



= MVT+Wf-l) -r 2 /2 



of the potential 



V(r) 



'^1 



The corresponding eigenvalue is 1 + + 87/2. We no- 
tice that the rate of decay increases as 7 becomes larger 
because the "intensity" of the positive singularity 7r -2 
also increases. 

If q > 2, then the decay at the singularity is exponen- 
tial. This is demonstrated by the decay at origin of the 



ground state eigenfunction 
1 



*o(r) 



■ exp 



of the potential 



V(r) 



1 



2 I r 2 ^ 1 ) 



7 



~ + "2 



7(7+1) 2 

r 2+7 r 7 



where g = 2(7 + 1). The ground state eigenvalue for the 
given example is 0.5. 

As the last example shows, the density matrix for a 
Lennard-Jones potential, potential that has a positive 
singularity of the form r -12 at origin, is decaying expo- 
nentially fast at r — 0. Because of this strong decay, the 
expectation values of the squares of the first and second 
derivatives of the potential against the diagonal density 
matrix are finite. As discussed in Section V.B, in these 
conditions we expect that the convergence rates of the 
reweighted techniques are cubic. A similar discussion 
holds for a one-dimensional system. For instance, the 
ground state eigenfunction of the potential 



V(x) = 



7(7+1) 



+£ 2 } , x > 0, 
x < 



is 



*o(x) = 6XP 



\x~> T 2 

0, 



x > 0, 
x < 



Again, ^o(x) has exponential decay at x = whenever 
7 > 0. (To recover the singularity of the Lennard-Jones 
potential, set 7 = 5). 



